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ABSTRACT: We study the predictability of large events in self-organizing 
systems. We focus on a set of models which have been studied as analogs of 
earthquake faults and fault systems, and apply methods based on techniques which 
are of current interest in seismology. In all cases we find detectable correlations 
between precursory smaller events and the large events we aim to forecast. We 
compare predictions based on different patterns of precursory events and find that 
for all of the models a new precursor based on the spatial distribution of activity 
outperforms more traditional measures based on temporal variations in the local 
activity. 
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I. INTRODUCTION 



Self-organized criticality (SOC) has received considerable attention over the past sev- 
eral years, as a possible means to explain scaling behaviors observed in a broad class of 
nonequilibrium systems including systems in geology, economics, and biology. 1 The theo- 
retical prototype is the sandpile model, in which sand is slowly added to a pile and released 
in instantaneous avalanches of a wide range of sizes which are triggered when the height 
(or slope or stress) locally exceeds a specified threshold. Self-organized criticality refers to 
the particular case of when the system size sets the cutoff for the largest events which are 
observed. More generally, self-organizing systems, whether critical or not, typically exhibit 
scaling over some range of sizes and are thought to evolve so that fluctuations in space and 
time are intrinsically coupled by an underlying threshold dynamics. There has recently 
been a considerable effort to use self-organizing systems as simple dynamical models of 
seismic phenomena, ''' in part due to the clear connection between earthquakes and 
threshold dynamics. Particular attention has been paid to the robust power-law scaling 
relation- the Gutenberg-Richter law 5 - relating the frequency of earthquakes to their size. 

An alternate direction of research concerns the predictability of the systems, which 
has important practical applications. The basic approach in such an endeavor is to utilize 
the available history of the system to forecast future events. Often one is most interested 
in predicting the largest events (e.g. great earthquakes), and it is the largest events with 
which we will be concerned here. Of course, in a well defined deterministic system precise 
knowledge about the present configuration of the system will yield very good, if not perfect, 
prediction. However, for many real systems specific equations describing the detailed 
evolutions of the systems are not known, and, in addition, the information on which the 
forecast must be based is typically incomplete. The basic prediction problem is, therefore, 
an inverse problem in the sense that one wants to use some information such as the time 
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series of events to infer something about the likely configuration of the system which is, in 
most practical situations, completely inaccessible to measurement. 

For example, earthquake catalogs list the date, time, location, and magnitude of de- 
tected events and thus provide one possible source from which one might hope to deduce 
information about local stresses on a fault. If correlations are detected they may lead to 
measurable precursors that are useful for forecasting. While certain seismicity patterns 
have been recorded in catalogs prior to some subset of the large earthquakes, in most cases 
the catalogs are too short to determine conclusively whether there is a statistically signifi- 
cant correlation between these patterns and large events. Dynamical models of earthquake 
faults can thus be particularly useful in the context of the prediction problem. Study of 
models allows us to consider catalogs of arbitrary size, from which we can make statisti- 
cally significant statements about both the intrinsic predictability of dynamical systems of 
this type and the value of current prediction algorithms. 

In this paper we address prediction issues in a variety of self-organizing systems. The 
algorithms that we use are similar to, and clearly motivated by, the work of Keilis-Borok, 
et a/. 6 which we describe briefly below. Our motivations for applying prediction algorithms 
to a broad class of systems is to try to ascertain what classes of precursory phenomena 
are consistently observed in all of the models. The point here is that no completely 
realistic dynamical model of faults presently exists, but if the real system resembles in 
any substantive way the threshold dynamics characteristic of the models, then precursors 
which are observed in a broad class of models may prove useful in real systems. In the 
course of this endeavor we develop a new type of precursor which is currently not in use 
in any form in the seismology community, which performs particularly well in the models. 
In addition, an unanticipated result from our simulations is that systems which have no 
apparent conservation law seem to be more predictable than those with a conservation law. 
Finally, while our findings may have important applications, we would like to point out 
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the modest nature of our results. The remaining open issues are extensive, with the most 
striking and presumably difficult issue being the development of an organized approach to 
deriving the most efficient precursors based on limited spatio-temporal information for a 
high dimensional dynamical system. In contrast, here we address the issue of predictability 
in the context of several fixed prescribed precursors used in isolation. 



II. THE MODELS 

We focus on a set of models which have been suggested as possible dynamical analogs 
of seismic phenomena which are described briefly below. We refer to the models as the the 
Bak, Tang, and Wiesenfeld (BTW) model, 1 the Olami, Feder, and Christensen (OFC) 
model, 2 the Chen, Bak, and Obukov (CBO) model, ^ and the Uniform Burridge and 
Knopoff (UBK) model. 4 ' 7 More detailed descriptions of the individual models may be 
found in the references cited above. 

We begin with the UBK model which satisfies a nonlinear wave equation 

d 2 U d 2 U 



U-<j>(U) + vt. (1) 



dt 2 dx 2 

Here U(x,t) represents the relative displacement of opposite sides of a homogeneous fault 
as a function of position x and time t. The variable v is the very slow shear rate driving 
the relative motion of the plates, and the key instability leading to chaotic behavior is a 
velocity-weakening, stick-slip friction law 4>(U). In the finite difference approximation, 
the model can be thought of as a one-dimensional chain of blocks, which is pulled slowly 
across a rough surface. 

Unlike the UBK model all of the other models ignore the details of inertial dynamics 
and friction laws, and instead evolve according to specified "breaking rules," so that when 
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the stress of a local block exceeds a threshold it relaxes according to some avalanche 
dynamics. For each of these the system can be thought of as a two dimensional 8 square 
lattice of "blocks" with open boundary conditions. In each case there is a particular rule 
which specifies the stress drop of the toppling site, the increases in stress of other sites, 
and the net stress drop of the system. 

The BTW model is the original sandpile cellular automaton. Of those we are consid- 
ering, it is the only model that is driven stochastically: on each iteration of the automaton 
the stress of a randomly selected site is increased by unity. If that site is above a specified 
threshold stress it initiates an avalanche in which each toppling site loses four units of 
stress, giving one unit to each neighbor (stress is dissipated at the boundary), and the 
cascade continues until all sites are below threshold. 

In the OFC and CBO models (as in the UBK model) stress is increased uniformly 
across the whole system. The OFC model is similar to the BTW model in that equal 
stress is transferred to each neighbor in each toppling, however, unlike the BTW model, in 
the OFC model the internal dynamics does not conserve stress. Instead, in each toppling 
the stress of the toppling site is set to zero, and each neighbor receives a fraction a < .25 of 
the initial stress of the toppling site (we will typically take a = .2). Any remaining stress is 
dissipated. This model has received considerable attention as a possible non-conservative 
example of SOC, although there has been some debate over whether this model exhibits 
SOC in the thermodynamic limit. ^ 

Like the BTW model, the CBO model does conserve stress away from the boundary. In 
addition, while it is not driven stochastically, it is does contain a stochastic element- after 
each site topples, its threshold stress is reset to a random value chosen uniformly from [0, 1]. 
In contrast, all of the other models have a fixed uniform threshold. Furthermore, among 
the models we are considering, the CBO model is unique in that it is the only one in which 
the relaxation dynamics explicitly takes place using long range elastic interactions, rather 
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than by redistributing stress only to neighboring sites. Here the stress of the toppling site 
is set to zero and the stress redistribution over the rest of the lattice is viewed as that due 
to a dipole force at the toppling site (hence falling off as r~ d ). 

While the BTW, OFC, and CBO models clearly differ from one-another in certain 
important ways, at least for the system sizes considered here they all generate pure power 
law event size distributions: 

P( s ) = (2) 

(see Fig. 1), analogous to the Gutenberg-Richter law 5 describing seismicity catalogs taken 
from the entire Earth or large regional fault systems. Thus we will refer to these systems 
as examples of SOC, with emphasis on "criticality" because the power law extends from 
the smallest event size up to essentially the system size. 

In contrast, the UBK model, while self-organizing, is not critical. The event size 
distribution consists of a power law describing the small to moderate events, and excess 
large events, which cut off at some characteristic size, independent of the system size 10 for 
systems which are large enough (Fig. 1). This statistical distribution is analogous to what 
is thought to apply to individual faults, or narrow fault zones, where the largest events 
appear to dominate the total slip, occurring at a rate which exceeds the extrapolated rate 
of small to moderate events. 11 ' 12 



III. THE PREDICTION ALGORITHM 



Our method of forecasting resembles the algorithm M8 introduced by Keilis-Borok 
et a/., 6 which is currently being studied as a possible means of using worldwide seismic 
data sets to predict the largest earthquakes in any given region. The M8 algorithm is 
based on the hypothesis that regional small scale seismicity may be used to diagnose an 
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upcoming large event. The procedure is to first coarse grain the catalogs in space and 
time, and then to measure certain precursors in these space-time windows. Several such 
precursors are used, and include a variety of measures based on the activity A, which 
within each space-time window is defined to be: A = ^earthquakes which are identified 
as main shocks 1 ^ and which are greater than or equal to some threshold size. Note that A 
is easily deduced from the time series of events in a region, and large values of A indicate a 
regional temporal clustering of events. More generally, an effective precursor is one which 
will typically sustain elevated (or depressed) values prior to a large event relative to its 
average value. 

In the Earth no single measure has yet been identified which reliably predicts all of 
the large events. Instead the M8 algorithm combines seven different precursors in a voting 
algorithm which is used to make predictions. It is important to articulate the prediction 
goal defined by Keilis-Borok et al. Instead of assigning some probability for an event to 
occur at a specific place and time in the future, the idea is simply to recognize certain 
seismicity patterns which might indicate a time of increased probability, or "TIP," for a 
large earthquake within a relatively large region in space (the spatial windows are chosen 
to be an order of magnitude larger than the event to be forecast). In particular, if a fixed 
number of precursors exceed individual thresholds in a region then the TIP is turned on. 

The earliest applications of these methods were based on existing data in real catalogs, 
and it was found that in order to capture roughly 80% of the large events, approximately 
20% of the total space-time volume had to sustain TIPs. Efforts are currently underway 
to evaluate the algorithms more thoroughly by establishing systematic tests for forward 
prediction. 14 However, use of the algorithm has been controversial for a variety of reasons 
including the intrinsic sensitivity to the inherent inaccuracies and incompleteness of the 
catalogs, 15 and sensitivity of the algorithm to features such as the initial placement of 
test regions (the space-time windows), where small adjustments in the spatial positions 
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of the regions and start dates of the catalogs can easily cause the algorithm to miss some 
of the events. 16 Of course, given a perfectly accurate catalog of arbitrary length (such 
as can easily be generated for models), the performance of these algorithms could easily 
be assessed. However, such catalogs are simply not available, so that the question of the 
predictability of earthquakes, as well as the development of effective algorithms, remain 
open and active areas of research. Here we will assess how well similar algorithms can be 
made to work on a collection of dynamical models. 

We consider a simplified version of this prediction algorithm, in which precursors are 
considered individually. We turn a TIP on when the precursor exceeds a threshold, and 
turn the TIP off when the precursor falls below the threshold. By varying the threshold 
we vary the total alarm time, and from this we construct a success curve 17 (Fig. 3), which 
plots the fraction of events predicted as a function of the fraction of the total space-time 
volume occupied by TIPs. 

Comparison of the success curves allows us to compare the effectiveness of different 
precursors as well as the predictability (based on these precursors) of different models. 
For a simple null test, our results can also be compared with the corresponding results 
for purely random methods, in which TIPs are issued completely arbitrarily. In that case, 
events are predicted purely by chance, and the fraction which are predicted successfully 
is simply given by the fraction of time the TIP is on: (% predicted) =(% alarm time), 
i.e. the success curve is the diagonal line. In a purely random system no algorithm will 
perform better than this method. In our case, this gives us an operational definition of 
predictability of models using a particular algorithm: if a model is predictable, the success 
curve for some precursor should deviate from this line in a statistically significant manner. 

It is not necessary for the success curve to lie above the diagonal line. In principal the 
curve could cross the line, or even lie completely below it. Any deviation from the diagonal 
is a sign that the precursor is detecting some correlation (or anticorrelation) in the system. 
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However, in practice to obtain a success rate which is better, rather than worse, than the 
corresponding results for random methods, one must base predictions on the complement 
of the original precursor whenever the success curve lies below the diagonal line. 18 

Finally, because of the limited amount of data which is available in seismic catalogs 
it is important to accompany any prediction with an assessment of the associated confi- 
dence level. In contrast, in our case we can generate catalogs that are arbitrarily long, 
and thus obtain results to an arbitrarily high level of precision. In particular, to verify 
that our success curves have converged to their asymptotic limits in time, we generate 
independent curves for a series of exponentially growing time intervals. This ability to 
check for systematic convergence is especially useful for models in which the predictability 
is marginal. 



IV. RESULTS 



We begin with the UBK model. A small segment of the catalog of events is illustrated 
in Fig. 2. For each event, a line is drawn through all the blocks which slip. While precursory 
small scale seismicity is clearly correlated with the epicenters of future large events (much 
more so for the model than for the earth) 19 it begins on average after half of the mean 
recurrence interval between large events has elapsed, so that from Fig. 2 alone it is not 
clear how accurate predictions based on this pattern will be. 

In Ref.[ 20] a detailed study of predictability of the UBK model revealed that among 
a set of precursors, the two most effective are activity A= # earthquakes, and a new 
measure which is a better measure of the development of spatial correlations. This new 
measure, which we call active zone size AZS, is defined to be: AZS = # blocks which 
have slipped (independent of the number of times) in each space-time window. With A 
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we are able to predict 90% of the large events, with alarms which occupy 15% of the 
total space-time volume (Fig. 3), which corresponds to alarms which occupy significantly 
smaller time intervals than the average duration of small scale seismicity prior to large 
events 20 in Fig. 2. However, the performance of AZS is even more impressive, leading 
to successful predictions of 90% of the large events when alarms occupy only 8% of the 
space-time volume. 21 In the UBK model the effectiveness of AZS can be traced to the fact 
that very little stress is relieved when a block slips in a small event. Instead small events 
serve as markers that the region is locally close to threshold. While the two precursors are 
clearly not independent, in contrast to A, AZS is a much more direct measure of the size 
of the region that is near the threshold for slipping, and thus ultimately leads to the more 
direct assessment of the probability of a large event. 

The question remains: which precursors which worked well for the UBK model are 
effective precursors for other self-organizing systems? To address this, we next consider 
the SOC models. Apart from some measurements of correlation functions between events 
of similar size in the OFC model 22 (which detected a tendency of large events to clus- 
ter in time) there has been little work to characterize the SOC models in terms of the 
predictability. 

As previously noted, the behavior of the UBK and SOC models differ from one another 
substantially. In principle, for the SOC models we could coarse grain the catalogs of 
events in space and time in a manner exactly like that used for the UBK model. However, 
because in the SOC models the largest events span essentially the entire system, and it is 
these largest events we wish to forecast, the most sensible choice is to define the spatial 
windows to correspond to the entire system. Furthermore, in these systems the distinction 
between small precursory events, and the large events which we attempt to predict is no 
longer a sharp feature, as apparent in the statistical distributions. For that reason, we 
set a somewhat arbitrary lower cutoff for events to predict which corresponds to the size 
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where we estimate (by eye) that finite size effects first become apparent (see Fig. 1). Our 
preliminary estimates of the predictability of small and medium size events indicate that 
in comparison the largest events are at least as predictable, and in most cases significantly 
more so than the others. 23 

The success curves for the SOC models, as well as our previous results for the UBK 
model, are illustrated in Fig. 3. For all of the models both A and AZS yield success curves 
which deviate systematically from the results obtained for random methods (the diagonal 
line), and thus lead to some measurable predictability. In each case along essentially the 
entire success curve AZS gives the greater deviation, and hence is more effective than A 
as a precursor for a coming large event. For the SOC models, we find that in most cases 
these measures are in fact anticorrelated with large events, with success curves falling 
below the diagonal. Whenever this is the case, we plot the complements of these measures, 
i.e. lack of activity A (quiescence) and lack of active zone size AZS, because these are the 
measures which would be used in practice to obtain a success rate which is better than 
random methods. When the complementary measure is used, TIPs are turned on when our 
previous measures A and AZS take values below, rather than above, a specified threshold. 
Of the SOC models, the OFC model is clearly most predictable, generating a success curve 
which is comparable to that of the UBK model. In this case the most effective measure 
is AZS, leading to 90% events predicted with alarm times of order 20%. Similarly, for 
the BTW and CBO models AZS outperforms measures based solely on activity (A for 
the BTW model, and A for the CBO model). However, compared to the UBK and OFC 
models, the gain over purely random methods is significantly reduced. 

In each case there is at least some correlation between small scale activity and coming 
large events, suggesting (but by no means proving) that self-organization may have impli- 
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cations for the predictability of real systems. The poor performance in the BTW and 
CBO models indicates that the correlations need not be strong, and may in real systems be 
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sufficiently weak that external noise or limited statistics may mask their presence. Recall 
that both the BTW and CBO models contain stochastic attributes, which we expect play 
important roles in limiting their predictability. In contrast, both the UBK and OFC mod- 
els, which exhibit the highest levels of predictability, are fully deterministic. However, the 
UBK and OFC models have an additional common feature which differentiates them from 
the BTW and CBO models. Both the UBK and OFC models do not satisfy a conservation 
law in the redistribution of total stress. A more detailed study is necessary to fully separate 
the relative roles of deterministic dynamics and the lack of conservation in determining 
the predictability of different systems, even within our limited definition of predictability. 
However, to address this question at least in part we have considered the predictability 
of the OFC model as the conservation parameter is varied. In that case we find that the 
predictability diminishes systematically as the level of conservation is increased, although, 
as illustrated in Fig. 4, even when the OFC model is fully conservative its predictability is 
clearly greater than that illustrated in Fig. 3 for the BTW and CBO models. 22,23 

In the SOC models the precursors based on quiescence are typically most effective 
because, unlike the UBK model, the stress on a site is set to zero each time a block slips, 
independent of the event size. Thus a lack of events is more likely to signify that the system 
is near the slipping threshold. Interestingly, in the Earth both increases and decreases of 
seismicity have been observed prior to large events, 25 suggesting that perhaps the UBK and 
SOC models may all contain some elements which are relevant to real faults. Ultimately 
it may be of interest to incorporate the more complete dynamical treatment of individual 
faults present in the UBK model into an SOC-type model which describes fault zones. 

In both the UBK and SOC models, large events involve large spatial regions, so that 
in order for a large event to occur, the system must be near threshold across a relatively 
large region in space. In all of the models considered here, methods based on AZS are 
more effective than those based on A, since AZS provides the more direct measure of 
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the development of such regions. The strong performance of this new precursor is par- 
ticularly noteworthy because such measures are not currently being used quantitatively 
in earthquake prediction algorithms such as M8. Certain tendencies towards clustering in 
space and time have been noted, 26 and are the phenomenological basis of these algorithms. 
However, the precursors which are used are based on measures such as A which within 
a given region track the development of temporal correlations. In such measures, spatial 
correlations are only accounted for in the most primitive way- i.e. in the initial definition 
of the spatial window. In the Earth including precursors which measure the development 
of geometric spatial correlations is complicated by the inhomogeneity of fault networks 
and the difficulties associated with accurately locating slip. Nonetheless, a box counting 
algorithm, in which the current spatial windows are coarse grained, and a count is made 
of regions exhibiting seismicity above or below some background level, may be adequate 
to measure the analog of AZS or AZS. It would be of significant interest to assess the 
performance of such a measure in comparison to activity based precursors in the Earth. 
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FIGURES 

(1) . Event size distributions P(s) vs. s for the (a) UBK, (b) OFC, (c) BTW, and (d) CBO 
models. In each case s is a measure of the size of the event- the integrated slip (seismic 
moment) for the UBK model, and the number of sites which topple for the others. In 
each case we attempt to predict events with s > s, where s is some characteristic size 
(in each case s is marked with an arrow). The UBK model exhibits a sharp distinction 
between small (s < s) and large (s > s) events, while the others exhibit power laws cut 
off at s ~ s (determined by finite size effects). In these cases, to exactly determine this 
crossover length a careful study of finite size effects must be performed. For the purposes 
of qualitative comparisons, it is sufficient to roughly estimate s as we do here. We take: 
N = 8192, a = .01, a = 3, and £/a = 10 for the UBK model, 20 system sizes 32 x 32 for 
the other models, and a = .2 for the OFC model. 2 

(2) . A small sample catalog as a function of space x (block number) and time t in the 
UBK model. A line segment marks blocks which slip in each event, and a cross marks the 
epicenter of each large event, which is clearly correlated with the small scale seismicity. 
The box corresponds to a space-time window within which A and AZS are evaluated. 

(3) . Success curves for the (a) UBK, (b) OFC, (c) BTW, and (d) CBO models. For each 
model, the best spatial measure (AZS or AZS) leads to more precise predictions than the 
best temporal measure (A or A). In each case, results for the complement measure (e.g. A 
vs. A) are obtained by a reflection of the curve across the diagonal. In the UBK model, we 
predict epicenters of large events, and optimal spatial windows are less than the size of the 
large event as shown in Fig. 2. We coarse grain the system into many overlapping spatial 
regions, and then obtain the success curve by calculating the fraction of events successfully 
predicted by some window containing the event vs. the total fraction of the space-time 
volume which is occupied bt TIPs. 21 In the other models, the spatial windows are taken 
to be the entire system, and the goal is to predict the large event. We have performed a 
crude optimization to select the time windows within which the precursors are evaluated. 
The values correspond to At = .1 for the UBK model (where the spatial windows were 
taken to be 213 blocks), 33 net grains added during the time window for the BTW model, 
.15 net stress added per site for the OFC model, and .007 net stress added per site in the 
measurement of AZS and .015 net stress added per site in the measurement of A for the 
CBO model. 

(4) . The deterministic and fully conservative OFC model. Here a = .25 whereas the results 
presented in Figs. 1 and 3 correspond to the dissipative case a = .2. Figure (a) illustrates 
the event size distribution, with s marked as in Fig. 1, and (b) illustrates the success 
curve, analogous to Fig. 3. Compared to the nonconservative case (Fig. lb), here the 
maximum event size is significantly increased, and compared to Fig. 3b the predictability 
is suppressed. As in Fig. 3, here we have crudely optimized over time windows to select a 
window which corresponds to .05 net stress added per site. 
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